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Abstract 

We study an individual-based predator-prey model of biological coevolution, using linear stability 
analysis and large-scale kinetic Monte Carlo simulations. The model exhibits approximate 1/f 
noise in diversity and population-size fluctuations, and it generates a sequence of quasi-steady 
communities in the form of simple food webs. These communities are quite resilient toward the 
loss of one or a few species, which is reflected in different power-law exponents for the durations of 
communities and the lifetimes of species. The exponent for the former is near —1, while the latter is 
close to —2. Statistical characteristics of the evolving communities, including degree (predator and 
prey) distributions and proportions of basal, intermediate, and top species, compare reasonably 
with data for real food webs. 
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I. INTRODUCTION 



Biological evolution presents many problems concerning interacting multi-entity systems 
far from equilibrium that are well suited for methods from nonequilibrium statistical physics 
[l|, 0] ■ Among these are questions concerningthe dynamics of the emergence and extinction of 
species on macroevolutionary timescales 0,0, 111 . Traditionally it has been common to treat 
ecological and evolutionary processes on very different timescales. However, it has recently 
been realized that evolution often can take place on short timescales, comparable to those 
of ecological processes jl, 0, 0, 0] • A well-known example of very rapid evolution is provided 
by the cichlid fishes of East Africa [Io|. Several models have therefore been proposed that, 
while spanning disparate scales of temporal and taxonomic resolution, consider the complex 
problem of coevolution of species in a fitness landscape that constantly changes with the 
composition of the community. Early contributions were simulations of parapatric and 



sympatric speciation [11] and the coupled NK model with population dynamics [1 
More recent work includes the Webworld model [§, 14, the tangled-nature model [1 
17 . 18| and simplified versions of the latter (l9l . I20I 21], as well as network models (22I . I23I 




Recently, large individual-based simulations have also been performed of parapatric and 



sympatric speciation [24|, |25| and of adaptive radiation 



Many of the models discussed above are deliberately quite simple, aiming to elucidate 
universal features that are largely independent of the finer details of the ecological interac- 
tions and the evolutionary mechanisms. Such features may include lifetime distributions for 
species and communities, as well as other aspects of extinction statistics, statistical prop- 
erties of fluctuations in diversity and population sizes, and the structure and dynamics of 
food webs that develop and change with time. 

In the present paper we continue our study of a simplified version of the tangled-nature 
model. In the early studies of this individual-based model of coevolutio n |19l . I20L l2l| , the 
interspecies interactions, which are described by an interaction matrix M [27J , were random 
and could produce any combination of pair interactions: favorable-favorable, deleterious- 
deleterious, or favorable-deleterious. Under those conditions the model was found to evolve 
through a sequence of quasi-stable communities, in which all species interact with mutu- 
ally favorable interactions, i.e., mutualistic or symbiotic communities. For that reason we 
shall hereafter refer to that version of the model as the mutualistic model. Here we instead 
concentrate on a version that specifically describes the evolution of predator-prey commu- 
nities. This restriction is enforced by means of an antisymmetric interaction matrix, so 
that an interaction that is favorable for one member of a pair is deleterious for the other, 
and vice versa. Many aspects of the dynamics of this predator-prey model are similar 
to the mutualistic model, such as approximate 1/f noise in species diversity and popula- 
tion sizes, and power-law distributions of the durations of communities and the lifetimes of 
species. However, some of the power-law exponents are different, and, most importantly, 
the predator-prey model produces communities that take the form of simple food webs. It 
also shares with the mutualistic model the property that mean population sizes and stability 
properties of fixed-point communities can be calculated exactly in the absence of mutations. 
Comparisons of some aspects of the predator-prey model (with a much smaller number of 
potential species than used here) to those of the mutualistic model have been presented in 



Refs. |28|, [29|. The focus of the present paper is a much more detailed discussion of the 
dynamics and the structure of the resulting food webs (including comparison with real food 
webs) for the predator-prey model with a large number of potential species. For this purpose 
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we use both exact linear-stability analysis and large-scale kinetic Monte Carlo simulations. 
In particular we wish to study the fluctuations in the statistically stationary state that de- 
velops for long times. A motivation is the hope that understanding of these stationary-state 
fluctuations can provide information about the system's sensitivity to external perturbations 
in a way analogous to a fluctuation-dissipation relation [30]. We therefore carry out very 
long simulations. 

The rest of this paper is organized as follows. The model is presented in Sec. [TTJ Ex- 
act linear-stability analysis is performed in Sec. IHH including fixed-point population sizes 
in Sec. IIII Al and stability considerations in Sec. IIIIBl Numerical results are presented in 
Sec. IIVt including species abundance distributions and time series of diversities and popu- 
lation sizes (Sec. IIV Aj) . power spectral densities (Sec. IIVBI) . species lifetimes (Sec. HVCj) . 



durations of evolutionarily quiet and active periods (Sec. HVDj) . and community structure 



and stability with comparison with real food webs (Sec. IIV Ep . Our conclusions are summa- 
rized in Sec. |V] and the method used to calculate the interaction matrices for systems with 
a large number of potential species is explained in Appendix |A] 

II. MODEL 

The model considered here is a version of the macroevolution model introduced by Rikvold 
and Zia 19[ as a simplification of the tangled-nature model of Jensen and coworkers [l6|, 



171 . |18||. In this version, the interspecies interactions are constrained to represent a pure 
predator-prey system. As in Ref. Ij|, selection is provided by the reproduction rates in an 
individual-based, simplified multispecies population-dynamics model with nonoverlapping 
generations. This interacting birth/death process is augmented to enable evolution of new 
species by a mutation mechanism. The mutations act on a haploid, binary "genome" of 
length L, as introduced by Eigen for molecular evolution [3l|, El. This bit string defines 



the species, which are identified by the integer label I G [0, 2 L — 1]. Typically, only a few of 
these 2 L potential species are resident in the community at any one time. 

During reproduction, an offspring individual may undergo a mutation that flips a ran- 
domly chosen gene (0 — > 1 or 1 — > 0) with a small probability, //. The mutation thus 
corresponds to diffusional moves from corner to corner along the edges of an L-dimensional 



hypercube |33|, 1341 ] . A mutated individual is assumed to belong to a different species than 
its parent, with different properties. Genotype and phenotype are thus in one-to-one corre- 
spondence in this model. This is clearly a highly idealized picture, and it is introduced to 
maximize the pool of different species available within the computational resources. This 
picture is justified by a large-scale computational study of the mutualistic version of the 
model studied in Ref. pjl], in which species that differ by as many as L/2 bits have cor- 
related properties [2lj]. Remarkably, this study reveals that the more realistic, correlated 
model has long-time dynamical properties very similar to the uncorrelated model. 

The reproduction probability Pi{t) for an individual of species I in generation t depends 
on the individual's ability to utilize the amount R of available external resources, and on 
its interactions with the population sizes nj(t) of all the species present in the community 
at that time. The dependence of Pi on the set of rij is determined by an interaction matrix 



M 27] with offdiagonal elements Mu that are continuously and symmetrically distributed 
on the interval [—1, +1] in a way defined specifically in the next paragraph. The elements of 
M are chosen randomly at the beginning of each simulation run and are subsequently kept 
constant throughout the run (quenched randomness). (For a discussion of how the matrix 
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elements are created for L > 13, in which case the 2 L x 2 L matrix does not fit into the 
memory of a standard workstation, see Appendix [A] This method leads to a distribution 
that is triangular on [—1, +1].) 

In contrast to our previously studied model (l9l . 20 . 21]. in which the interaction matrix 



has no particular structure, the predator-prey dynamics is enforced by the requirement that 
the off-diagonal part of M must be antisymmetric. Thus, if Mjj > and Mjj < 0, then 
species I is the predator and J the prey, and vice versa. In order to keep the connectance 
of the resulting communities consistent with food webs observed in nature 35, 36[ the 



(M/j, Mjj) pairs are chosen nonzero with probability c = 0.1. The nonzero elements in the 
upper triangle of M are chosen independently from the triangular distribution on [—1, +1], 
described in Appendix |A] Self-competition is included in the model by choosing the diagonal 
elements of M randomly and uniformly from [—1,0). 

The reproduction probability for species /, Pi(t), depends on R and the set {nj(t)} 
through the nonlinear form, 

Pl{t) = l + exp[-A,(JUM*)})] ' (1) 

where 

Aj(R, {nj(t)}) = -h + Vl R/N tot + Mjjnj{t)/N tot . (2) 

j 

Here bj is the "cost" of reproduction for species I (always positive), and rjj (positive for 
primary producers or autotrophs, and zero for consumers or heterotrophs) is the ability of 
individuals of species I to utilize the external resource R. The latter is renewed at the same 
level every generation and does not have independent dynamics. The total population size 
is N tot (t) = J2j ti j(J : )- P 11 contrast, the total number of species present in generation t (the 
species richness) will be defined as A/"(i).] The population-limiting reproduction costs bj are 
chosen randomly and uniformly from the interval (0, +1]. Only a proportion p of the 2 L 
potential species are producers that can directly utilize the resource. (For the numerical 
data shown here, we use p = 0.05.) Thus, with probability (1 — p) the resource coupling 
rjj = 0, representing consumers, while with probability p the rjj are independently and 
uniformly distributed on (0, +1], representing producers of varying efficiency In addition 
to the constraints on M mentioned above, we require that producers always are the prey of 
consumers. Thus, the case rjj > and rjj = with Mjj = —Mjj > is forbidden and is 
changed during setup of the matrix by reversing the signs of the interactions for the pair in 
question. 

For large positive A/, (small birth cost, strong coupling to the external resources, and 
more prey than predators), the individual almost certainly reproduces, giving rise to F off- 
spring. In the opposite limit of large negative A/, (large birth cost, weak or no coupling to 
the external resources, and/or more predators than prey), it almost certainly dies without 
offspring. The nonlinear dependence of Pi on A/ thus limits the growth rate of the popula- 
tion size, even under extremely favorable conditions. It also sets a practical negative limit 
on Aj, below which conditions are so unfavorable that reproduction is virtually impossible. 
(A more general version of Eq. (j2J), in which population growth is directly limited by a "Ver- 



hulst factor" [37[ or "environmental carrying capacity" [381 ] as is necessary in models that 
allow mutualistic interactions, is discussed in Ref. [29(.) We note that "energy dissipation" 
in this model is achieved through the birth costs bi and the self-competition terms Mjj. 
Equivalent effects could have been produced by making the positive Mjj smaller than the 
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corresponding negative ones by an "ecological efficiency" factor between zero and unity, as 



in the Webworld model 18, 14, 15 



The normalization of A/ with N to t implies global competition. This is not very realistic, 
but it enables us to find exact expressions for the stationary values of the average population 
sizes in the mutation-free limit. (See Sec. lIII Al ) The model can thus be used as a benchmark 
for more realistic ones in future research. 

An analytic approximation describing the development in time of the mean population 
sizes (averaged over independent realizations), (rij{t)), can be written as a set of coupled 
difference equations, 

(n 7 (t+l)> = (n/(t)>FP / (i2,{(n J (t)>})[l-^] 

+(l,/L)Fj2(nK(i)(t))PK(i)(R, {(nj(t))}) , (3) 

K(I) 

where K(I) is the set of species that can be generated from species / by a single mutation 
("nearest neighbors" of / in genotype space). 

III. LINEAR STABILITY ANALYSIS 
A. Fixed-point communities 

An advantage of the model studied here is that its fixed-point communities in the 
mutation-free limit can be found exactly within a mean-field approximation based on Eq. ([3]) 
(l9| . To obtain a stationary solution for a community of M species, we must require Pj = 1/F 
for all M species. Equations ([I]) and (T5]) then give rise to M linear relations, which can be 
written on the matrix form 

-\b)N: ot + \ v )R + M\n*)=0, (4) 

where bj = bj — ln(P — 1), \b), \r)), and |n*) are the column vectors of bi, rji, and n*j, 
respectively (in all cases including only those M species that have nonzero n*j), and M is 
the corresponding M x M submatrix of M. (For simplicity, we drop the ( ) notation for the 
average population sizes, and the asterisk superscripts denote fixed-point solutions.) 
The solution for \n*) is 

\n*) = -M- 1 [\ V )R-\b)N: ot ] , (5) 

where M _1 is the inverse of M. To find each n}, we must first obtain A^ t * ot = (1|^*), where 
(1| is an A^-dimensional row vector composed entirely of ones. Multiplying Eq. ([5]) from the 
left by (1|, we obtain 

BE + 9iV t * ot = , (6) 



where the coefficients 



e= i-<ijft-|*> and £= (Hft-'h>> (7) 

(1IM- 1 !!) (IIM- 1 )!) 
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have been written with (1|M _1 |1) in the denominators in order to remain finite even for 
near-singular M [2^] . They can be viewed as an effective interaction strength and an effective 
coupling to the external resource, respectively. The solution of Eq. (jSj) is 

__R£ R(1\M~'\ V ) 
tot " 6 - (l|M-!|fe)-r W 

To find each n*j separately, we now only need to insert this solution for jV t * ot in Eq. (J5j). 

Only those |n*) that have all positive elements can represent a feasible community 39 . 
If M = or is otherwise singular, the set of equations (jl]) is inconsistent for M > 1, unless 
bj and 77/ both are independent of / (this case is equivalent to Af = 1). The only possible 
stationary community then consists of one single species, the one with the largest value of 



r\ijbj. This is a trivial example of competitive exclusion [40j, |4l|, |42j. If rji/bi has the same 
value for all M values of J, we have an example of a neutral model 43 . 

In Ref. 2^1 it was shown that Eq. ([6]) for fixed £ and G can be seen as a maximization 
condition for a "community fitness" function, 



(9) 



This result would not be particularly remarkable if £ and G were externally fixed parameters. 
However, extinctions and mutations provide a mechanism for both parameters to change as 
old species go extinct and new species emerge. In our numerical simulations we find that 
their values evolve toward and then fluctuate around values that maximize $, limited only 
by the internal constraints on M and \bj). In particular, this means that G approaches 
closely to from the negative side [29]]. In the simulations presented in this paper (see 
Sec. IIVI) averages over the final communities of twelve independent simulation runs yield 



G « —0.10 ± 0.03 and £ w 0.7 ± 0.2, corresponding to N* ot « (7.4 ± 1A)R. In comparison, 
for fourteen random, feasible communities obtained as described in Sec. IIVE II below, the 
corresponding parameters are Q ~ —0.39 ± 0.06 and £ « 0.66 ± 0.08, corresponding to 
N* ot « (2.0 ± 0.3)R. A detailed discussion of the statistical properties of these and other 
quantities characteristic of the simulated communities are given in Sec. IIVE 21 



B. Stability of fixed-point communities 



The internal stability of an A/"-species fixed-point community is obtained from the matrix 
of partial derivatives, 

dnj{t + l) 



dnj(t) 



Su + A 



1.) 



(10) 



|n*> 



where 5u is the Kronecker delta function and A j, j are elements of the community matrix A 

R Vl + (M\n*})j 



381 ] . Straightforward differentiation yields [29j 



A 



1. 1 



1 

1 

F 



nj 

AT* 
iv tot 



M u - 



N* 

iv tot 



(11) 



where {jS/l\n*))i is the element of the column vector M|n*}, corresponding to species /. In 
order for deviations from the fixed point to decay monotonically in magnitude, the magni- 
tudes of the eigenvalues of the matrix of partial derivatives in Eq. ffTUl) . A + 1, where 1 is 
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the A/"-dimensional unit matrix, must be less than unity. The value of the fecundity used in 
this work, F = 2, was chosen to satisfy this requirement for J\f = 1. 

Since new species are created by mutations, we must also study the stability of the fixed- 
point community toward "invaders." Consider a mutant invader i. Then its multiplication 
rate, in the limit that <C nj for all TV species J in the resident community, is given by 

nAt+l) F . . 

— = . (12) 

m{t) 1 + exp [-\{R, {n}})} 1 ; 

The Lyapunov exponent, In Wij ( t + 1) / Uj (t)], is the invasion fitness of the mutant with respect 
to the resident community 44, |45[. It will be studied numerically in Sec. HVEll 



IV. NUMERICAL RESULTS 

We performed twelve independent, long simulation runs of 2 25 = 33 554432 generations 
of the model with the following parameters: genome length L = 20 (1048 576 potential 
species), external resource R = 2000, and mutation rate p, = 10~ 3 , with connectance pa- 
rameter c = 0.1 and a proportion p = 0.05 of the potential species as producers. These 
parameters were chosen to represent the realistic situation that the number of species res- 
ident in the community at any time is much smaller than the number of potential species 
(i.e., that Af(t) <C 2 L ), and also that Af(t) <C N tot (t) so that at least one species has a sub- 
stantial population size. In this parameter range the model is not very sensitive to the exact 
parameter values j2^|. We also note that the chosen connectance is above the percolation 



limit on the L = 20 dimensional cube of potential genotypes [34J, |46[], so that a finite fraction 
of the genotypes can be connected by mutations along paths of nonzero interactions. 

The very long simulation times were chosen because our main interest is in the statistically 
stationary dynamics of macroevolution over timescales much longer than the ecological ones 
of a few generations. Each run therefore starts with a "warm-up period" of about one 
million generations before the 2 25 -generation data-taking period. Details of the simulation 
algorithm were given in Ref. [HI. 



A. Time series and species abundance distributions 

We collected time series of a number of quantities including several measures of diversity 
or species richness, as well as population sizes of producer and consumer species. Time series 
of diversities and population sizes for one representative run are shown in Fig. [IJ 

In order to filter out noise from the small-population species that are mostly unsuccessful 
mutations, we use the diversity measure known in ecology as the exponential Shannon- 
Wiener index [47j . It is defined as the exponential function of the information-theoretical 
entropy of the population distributions, D(t) = exp [S ({n>i(t)})], where 

S({nj(t)}) = - Yl Pi(t)hPi<t) ( 13 ) 

{/| PI (i)>0} 

with pi(t) = ni(t)/N tot (t) for the case of the curves labeled "All species" in Fig. [lj For the 
producers or consumers, the sums and normalization constant include only the appropriate 
species. The utility of the Shannon- Wiener index is illustrated by the data presented in 
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FIG. 1: (Color online.) Time series for the Shannon- Wiener diversity index (a) and the population 
sizes (b). The jagged curves in the background show data sampled every 2048 generations to 
show the rapid fluctuations, while the curves in contrasting color/brightness in the foreground 
are running averages over 524 288 generations, emphasizing the slower fluctuations. The three 
sets of curves represent all species (black/light gray; black/yellow online), producer species (light 
gray/dark gray; green/violet online), and consumer species (medium gray/light gray; red/cyan 
online). 



Fig. El In Fig. 12(a) we show two versions of the Species Abundance Distribution (SAD) for 
the model. This is one of the tools most widely used in ecology to describe the distribution 
of the number of species over population sizes j48|. The SADs shown as full lines in the 
figure represent full communities, from which we have only removed any consumer species 
that are not connected to the external resource through an unbroken chain of nonzero 
interactions. These we term "full connected communities" ("full communities" for short). 
(The removal of disconnected consumer species actually has a numerically insignificant effect 
on the SADs.) The SADs for the full communities are dominated by a large number of species 
with very small populations, which to a large extent represent unsuccessful mutants. This 
was explicitly shown by extracting "core communities" in the following way. Communities 
were sampled every 256 generations, and species with population sizes below eight were 
excluded. It was then checked whether each included species also existed with at least 
this minimum population 256 generations ago, and if this was not the case, the species was 
removed from the community as unstable. The fixed-point populations for the community of 
remaining species were then calculated according to Eq. (jSJ), species with negative fixed-point 
populations (unfeasible species) were removed, and the fixed-point calculation was repeated 
until all species remaining in the community had positive populations. The SAD was then 
calculated for each of these feasible core communities and averaged over all communities in 
the twelve independent simulation runs. This procedure removes most of the low-population 
species, as shown by the dashed curves in Fig. 121(a). The core-community SAD appears to 
be intermediate between Fisher's log-series distribution jH, El and Preston's log-normal 



pi, and (3 in the function [51 



distribution [48|, |5fJ]. It can be semiquantitatively approximated by fitting the constants C, 

p(n) = C^n^e-^/TiP) , (14) 
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FIG. 2: (Color online.) (a) Species abundance distributions (SADs), normalized to the number 
of species. Solid curves represent "full connected communities," and dashed curves represent "core 
communities," both extracted as described in the text. Both were sampled every 256 generations. 
The data were averaged over twelve independent simulation runs, and the error bars represent 
standard errors, based on the differences between runs. The dot-dashed curve is a fit to the curve 
describing all core species, using Eq. (jTSj) with parameters C = 13.569, fi = 0.00210967, and 
(3 = 1.89338, interpolating between log-series and log-normal forms, (b) Histograms of the full- 
community species richness, the species richness of the core communities, and the Shannon- Wiener 
diversity for the sampled communities. The latter is seen to be an excellent approximation for the 
species richness of the core communities. 



which interpolates between these two limiting forms. The fit is shown by the dash-dotted 
curve in the figure. Additional evidence for the agreement between the Shannon- Wiener 
diversity index and the species richness of the core communities is shown in Fig. [2](b), where 
histograms of the two are in excellent agreement and both show a narrow peak near fifteen 
species (mean values of 15.7 and 15.3 species, respectively), while the raw species richness 
yields a wide distribution with a mean of 49.3 species. All three distributions are very 
well fit by gaussians and are much more symmetric than diversity distributions arising from 



Rossberg et al.'s speciation model of food webs [52 



Time series of additional quantities that indicate the level of evolutionary activity are 
shown in Fig. [3] for the same simulation run as in Fig. [TJ These are the magnitude of the 
logarithmic derivative of the Shannon- Wiener diversity index (Fig. [3](a)), and the size of 
extinctions per generation (Fig. [3](b)). The latter is denned as the sum of the maximum 
population sizes reached by the species that go extinct in each generation. 

The properties of the fluctuations in the time series were analyzed with several methods, 
including power spectral densities (PSD), lifetimes of individual species, and the durations 
of quiet and active periods during the evolution. The results of each are reported in the 
following subsections. 



B. Power spectral densities 



PSDs of the diversity and population-size fluctuations, averaged over the twelve indepen- 
dent simulation runs, are shown in Fig. HI both for the total population and for the producers 
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FIG. 3: (Color online.) Time series for the same simulation run shown in Fig. [H displaying 
quantities that measure the evolutionary activity: the magnitude of the logarithmic derivative of 
the Shannon- Wiener diversity index, \dS/dt\ (averaged over 16 generations) (a), and the size of 
extinctions per generation (b). The horizontal lines in each part indicate cutoff levels used to 
define quiet and active periods as discussed in Sec. HVDi 

and consumers separately. The spectra indicate 1/ / like noise over more than five decades in 
time. A weighted fit to the PSD for the overall diversity in Fig. [U^a) yields a power law f~ a 
with a m 1.29 ± 0.01. This power is also seen to fit reasonably well, both with the data for 
the diversities of producers and consumers over the whole frequency range in Fig. H^a), and 
with the PSDs of all three population measures at low frequencies in Fig. IU(b), as well as 
for the extinction measures at low frequencies in Fig. IH(c). This suggests that the long-time 
fluctuations in the diversity, as well as in the population sizes and the extinction measures, 
obey the same power law on long timescales. On short timescales the PSDs for the popula- 
tion sizes have a more complicated structure, possibly indicating overdamped oscillations on 
a scale of a few hundred generations. The extinction measures show a wide region of white 
noise for high frequencies, due to the frequent extinction of unsuccessful mutants. However, 
the behaviors for low frequencies appear consistent with the diversities and the population 
sizes. 



C. Species lifetimes 



The statistics of the lifetimes of individual species are characteristic of the evolution 
process. The species lifetime is defined as the time from a particular species enters the 
community, till it goes extinct (i.e., its first return time to zero population size). Histograms 
showing the distributions of species lifetimes, for all species as well as for producers and 
consumers separately, are shown in Fig. Although there are some undulations in these 
curves, they remain close to a power law t~ T1 with exponent T\ ~ 2 over more than six 
decades in time, which is the maximum we could expect with the length of our simulations. 

The t~ 2 dependence of the lifetime distributions is quite universal. It is found, e.g., in 



our previous studies of the mutualistic model [2ll . l28l . |29J , and it is in general characteristic 
of stochastic branching processes [53j |. 
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FIG. 4: (Color online.) PSDs for the Shannon- Wiener diversities (a), population sizes (b), and 
the extinction sizes and number of species going extinct per generation (c), all averaged over 
twelve independent runs of 2 25 generations. The error bars are standard errors, estimated from 
the variations between the individual runs. The dashed straight line with slope —1.29 in (a) is a 
weighted fit to the PSD for the overall diversity over the whole frequency range. The dashed lines 
in (b) and (c) are guides to the eye with the same slope. 

Lifetime distributions for marine genera that are compatible with a power-law exponent 
in the range —1.5 to —2.5 have been obtained from the fossil record 0,13, H|- However, the 
possible power- law behavior in the fossil record is only observed over about one decade in 
time - between 10 and 100 million years - and other fitting functions, such as exponential or 
log-normal, are also possible. Nevertheless, it is reasonable to conclude that the numerical 
results obtained from complex, interacting evolution models that extend over a large range 
of time scales support interpretations of the fossil lifetime evidence in terms of nontrivial 
power laws. 
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FIG. 5: (Color online.) Histograms for the lifetimes of individual species. The dashed, straight 
line is a guide to the eye, corresponding to a t~ 2 power law. Data averaged over 12 independent 
simulation runs. 




FIG. 6: (a) Histogram of the logarithmic derivative of the overall diversity, dS/dt. The data were 
averaged over 16 generations for each run and then over 12 independent runs. The dashed curve 
is a gaussian fit. The time series of the magnitude of this quantity for one particular simulation 
run is shown in Fig. [3ja). (b) Histograms for the durations of active and quiet periods for various 
values of the cutoff y c for \dS/dt\, averaged over 12 independent runs. The dashed, straight line 
with slope —1.07 is a weighted fit to the histogram for the quiet periods for y c = 0.010 between 10 
and 10 6 generations. The steep histogram curve with only three data points between 10 and 100 
generations corresponds to the active periods, which are always very short. 



D. Quiet and active periods 



From the time series shown in Figs. Q] and [3] one sees that periods of moderate fluctu- 
ations are punctuated by periods of high activity. The communities corresponding to the 
lower fluctuation intensities are known as quasi-steady states (QSS). (In the literature on 
the tangled-nature model 0, E3 the QSS are referred to as quasi evolutionarily steady 
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FIG. 7: Histograms of the durations of quiet periods, obtained from the time series of extinction 
sizes, an example of which is shown in Fig. Mb). Data averaged over 12 independent simulation 
runs. The straight, dashed line with slope —1.07 is a guide to the eye, based on the fit to the QSS 
duration distribution based on \dS/dt\, shown in Fig. 0(b). 

strategies, or q-ESS.) The cores of these communities correspond to the fixed-point com- 
munities of the mutation-free system [lj|, 20], as discussed in Sec. IIV Al One measure of 



the degree of activity is the time derivative of the entropy or, equivalently, the logarithmic 
derivative of the total Shannon- Wiener diversity. Its magnitude is shown as a time series 
in Fig. [3](a), and a histogram is shown in Fig. [6](a). While the central part of the distri- 
bution is well approximated by a gaussian, the heavier wings appear to be exponential or 
even power-law. Quiet and active periods were defined as contiguous periods during which 
\dS/dt\ stayed below or above a cutoff y c , respectively. Octave-binned histograms for the 
probability distributions of the durations of quiet and active periods are shown in Fig. Mb) 
for various cutoffs. For the quiet periods a power law is seen with exponent near —1 (a 
weighted fit between 10 and 10 6 generations gives t~ T with r « 1.07 ±0.01) and a long-time 
cutoff that increases with increasing y c . The active periods for all values of y c are very brief 
in comparison. (Their histogram for y c = 0.010 is the steep curve with only three data 
points between 10 and 100 generations in Fig. E(b).) As a consequence, the system spends 
most of its time in QSS communities - a situation consistent on the community level with 



Eldredge and Gould's concept of punctuated equilibria [54], |55j, |56| . 

Durations of quiet periods could also be obtained from the time series of extinction sizes 
in Fig. Mb). Due to the white noise at short timescales, which was also apparent in the 
PSDs in Fig. Sic), the power-law behavior is limited to a window of longer times between 
about 1000 generations and the strongly cutoff-dependent long-time decay. As a result, this 
quantity does not provide as clear a quiet-period distribution as the entropy derivative. We 
therefore did not perform any independent fit to obtain a power-law exponent for the QSS 
durations measured this way. See Fig. [71 The decay with time is qualitatively consistent 
with that observed in Fig. Mb) for the QSS duration distributions based on \dS/dt\. 

It is quite remarkable that the exponent for the QSS durations is significantly different 
from the one for the species lifetimes, which is close to 2. This is particularly so because 
the two exponents appear to be approximately the same (both near 2) for the mutualistic 
version of the model [21, 28, 2^|. We believe that the explanation lies in the structure of the 
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FIG. 8: (Color online.) Major populated species shown vs time for the same simulation run 
shown in Figs. [T] and [3j The horizontal lines correspond to the species label, and the grayscale 
(color online) to the population size. Black: ni > 1000. Gray (red online): n/ € [101,1000]. (a): 
Producers, (b): Consumers. 



QSS communities generated by the present evolution process, which take the form of simple 
food webs. These are studied in Sec. IIVEI below. 



E. QSS community structure and stability 

1. General considerations 

The evolution process in the present model generates dynamic communities, in which 
species emerge, exist for a shorter or longer time, and eventually go extinct. The emergence 
and extinction of a major species are quite fast processes on the evolutionary timescale, 
and so the vast majority of randomly selected communities are QSS communities. This is 
confirmed by the short durations of evolutionarily active periods, shown by the corresponding 
histogram in Fig. [6](b). Diagrams of the population sizes of major producer and consumer 
species as functions of time in a particular simulation run are shown in Figs. E^a) and 
(b), respectively. In these figures a horizontal line represents a species. The beginning 
of the line represents the emergence of the species, and the end represents its extinction. 
The population size is represented by the color. We see that some species persist for tens 
of millions of generations, while others are so short-lived as to hardly be visible on the 
scale of these figures. This is consistent with the power-law behavior of the species-lifetime 
distribution (see Fig. E]). We also see that producer species appear to emerge and go extinct 
relatively independently of each other, while there is a significant correlation between the 
producer and consumer species. This correlation indicates that extinction of a producer 
species is likely to trigger a (limited) cascade of consumer extinctions. Conversely, a new 
producer species is likely to quickly acquire a group of consumer species. The structure of 
the plots in Fig. [8] contrasts with that of similar plots for the mutualistic model of Ref. 1^ 



(see Fig. 2 of that paper), in which species tend to emerge as well as go extinct together. We 
believe this is the reason for the difference between the exponents for the species lifetimes 
and the durations of QSS communities in this model: the overall community is relatively 
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FIG. 9: (Color online.) Food webs representing QSS core communities for the simulation shown 
in Figs. Q3 El and [8] at times near 22 x 10 6 generations (a), near 27 x 10 6 generations (b), and at 
the end of the simulation near 35 x 10 6 generations (c). The core communities were identified as 
described in Sec. IIV Al The thickness and head size of the arrows correspond to the magnitude of 
Mjj, and the area of the circles to the stationary population size, as calculated analytically from 
Eq. (0). Light gray (red online) lines connect nearest neighbors in genotype space. Labels in red 
mark species that persist from the first to the last snapshot, (a) to (c), while blue labels mark 
species that persist from (a) to (b), and and green labels mark species that persist from (b) to 
(c). Black labels indicate species that are unique to the particular community. 



resilient toward losing or gaining a single species [g, 35, 57, 58[. As a community it is more 
long-lived than the individual species, leading to the smaller value for the exponent of the 
distribution for the QSS durations (compare Figs. [5] and E](b)). (A similar relationship has 
been noted between the lifetimes of orders and their constituent genera in the fossil record 

Three representative QSS core communities from the same simulation run shown in 
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FIG. 10: (Color online.) Histograms of the multiplication ratio, rii(t + l)/nj(t) (exponential of the 
invasion fitness, Eq. (|12p) for low-density invaders. Each curve is averaged over 14 different core 
communities. Black: nearest neighbors in genotype space against stable QSS communities. Dark 
gray (red online): All species not in the community against stable QSS communities. Light gray 
(orange online): All species not in the community against random feasible communities. Inset: 
The part of the distributions for potentially successful invaders, rii(t + 1) / rii(t) > 1, with the y-axis 
on logarithmic scale. The tails appear nearly exponential. 



Figs. (U [3j and [8] are shown in Fig. [9j These are consecutive core communities near 22 x 10 6 , 
27 x 10 6 , and 34 x 10 6 generations, respectively, and they are all stable by the eigenvalue 
criterion discussed in Sec. IIII Bl The communities take the form of simple food webs with 
two trophic levels above the resource node. Communities with three trophic levels are also 
occasionally observed. The different branches of the webs are relatively independent of each 
other, and many consumer species have more than one prey. (See further quantitative dis- 
cussion in Sec. IIVE2I ) Both features contribute to the resilience against mass extinctions 
discussed above. 

The relative stability of evolved QSS core communities with respect to invasion by mu- 
tants is illustrated in Fig. [TUJ The results are averaged over 14 QSS core communities - 
the three shown in Fig. [9] plus the final communities of the eleven other simulation runs. 
This figure shows the multiplication ratio of a small population of invaders (the exponential 
function of the invasion fitness), given by Eq. f[T2"j) . Only about 2.3% of species outside the 
communities have multiplication ratios larger than unity, and most of these lie between 1.0 
and 1.1. This percentage does not seem to depend significantly on the Hamming distance 
of the invader from the community. We note that the species that are removed from the 
community during extraction of the core community are among these low invasion fitness 
species. This is a further indication that the difference between core and full communities 
is mostly made up by unsuccessful mutants. 

We also tested if the evolved core communities are more resilient toward invasion than 
randomly constructed feasible communities. Such communities are more difficult to con- 



struct in this model, than in the mutualistic model [19|. To have the same level of statistics 



as for the QSS communities, we produced 14 such communities in the following way. We 
started a run with a random sample of 200 species, each with ni = 10, and evolved the 
community for 1024 generations without mutations. We then tested the remaining commu- 
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FIG. 11: (Color online.) Log-linear plot of histograms of the birth cost bj (left) and producer 
coupling to the external resource, r\i (right), for members of QSS core communities (shaded areas) 
and full communities (lines only) . Both parameters are selected away from the uniform distributions 
of the full species pool (represented by the horizontal dashed line). 



nity for feasibility and removed species with a negative stationary population according to 
Eq. (jSJ). The resulting feasible communities had much smaller diversities than the evolved 
QSS communities - an average of only 3.6 species per community. These communities had 
a significantly larger percentage of potential invaders - about 4%, with a maximum multi- 
plication ratio near 1.7 (see inset in Fig. [TOi) . While there thus is a clear difference between 
the stability against invaders of QSS communities and random feasible communities, the 
difference is not as large as it is for the mutualistic model of Ref. (See Fig. 3 of that 
paper.) 



2. Detailed community structure and comparison with real food webs 

Next we turn to a detailed statistical description of the QSS core and full connected com- 
munities identified in Sec. IIV Al Since communities were extracted every 256 generations 
from twelve independent runs of 2 25 generations each, data were averaged over 131 072 x 12 
communities of each type. However, due to the long-time correlations in the evolution 
process, many of the communities from the same run are identical or similar. This sam- 
pling method thus ensures that the statistics for both communities and individual species 
are weighted according to their longevities. First we consider the properties of individual 
species, and next we turn to the collective properties of the corresponding food webs and a 
comparison with data for real food webs. 

The individual species are characterized by the birth cost bi, the self interaction Ma, and 
for producers by the resource coupling t]i. In the total pool of potential species these are 
uniformly distributed on (0, +1], [—1, 0), and (0, +1], respectively. Not surprisingly, the most 
prevalent species in the core communities turn out to be the most "individually fit" ones 
in the sense that they have low birth cost and, for producers, relatively strong coupling to 
the external resource. The resulting probability densities for bi and r\i are shown in Fig. [TTJ 
The selection for low birth cost is very strong for members of the core communities, and less 
so when the full communities are considered. This is further indication that the species that 
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FIG. 12: (Color online.) Histograms of community-level parameters, (a) The interspecies inter- 
action strengths, |M/j|, for QSS core and full communities, (b) The effective interaction strength 
(left) and resource coupling £ (right) for QSS core communities, defined in Eq. (JZJ). 

are ignored when extracting the core communities have higher bi, and thus overall are less 
individually fit, than the core species. On the other hand, the moderate selection for strong 
coupling to the external resource (large rjj) is approximately equal for producers in both 
types of communities. In contrast, the self interactions appear to be evolutionarily neutral, 
and their probability density remains approximately uniform for the members of both core 
and full communities (not shown). 

The first community-level quantities we consider are the nonzero interspecies interaction 
strengths, |Mjj|. Histograms for these, based on the same sampling and averaging as previ- 
ous results, are shown in Fig. [T27 a). The distribution for the core communities is significantly 
skewed toward strong interactions, compared to the triangular distribution characteristic of 
the total species pool, which is shown as a dashed line. The bias is much weaker when all 
members of the full communities are considered. (See Appendix [S] for a discussion of the 
unbiased distribution.) 

The biased distribution of interaction strengths combines with the strongly skewed dis- 
tributions of bj and rjj to produce values of the effective interaction strength and effective 
resource coupling £ for the core communities (defined in Eq. (j7j)) that are in excellent agree- 
ment with the expectations expressed in Sec. IIII Al As shown in Fig. [TWb). G is narrowly 
distributed closely below (6 = —0.10), while £ has a broader distribution over positive 
values with mean £ = 0.71. These averages are in excellent agreement with those obtained 
from the final communities of the twelve runs, obtained in Sec. IIV A[ 

Next we consider several quantities from network theory that can be compared with 
data from real food webs. These include the directed connectance, C = L/S 2 , where L is 
the number of links (i.e., the number of nonzero Mjj, Mjj pairs) and S is the diversity 
(here defined as the species richness), the linkage density, z = L/ S = CS, the probability 
distributions of a species' number of prey species (its generality or indegree), number of 
predator species (its vulnerability or outdegree), and their sum (its total degree) (52| . 

For comparison with the simulated food webs generated by our evolutionary model, we 
use data for seventeen empirical webs, including both aquatic and terrestrial communities. 
These data sets were kindly provided to us by J. A. Dunne. Tabular overviews of the 
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FIG. 13: (Color online.) Histograms for the linkage density, z = L/S (main part) and the directed 
connectance, C = L/S 2 (inset), for the model communities. 



parameters characterizing these communities are available in the literature |35|, l57|, [59|, [60 . 
For all the empirical webs we defined the diversity as the species richness S in terms of trophic 
species. These are obtained by lumping together all taxa that share all their predators and 



prey |59|, |61]. Like the simulated communities analyzed, the empirical webs contain no 
disconnected species or subwebs |59l| . The included empirical food webs are: Coachella 
Valley [62|, El Verde Rainforest [63l |. Scotch Broom [HII], St. Martin Island [HHI, and U.K. 
Grassland [HH1 (terrestrial); Bridge Brook Lake [13], Little Rock Lake [HI], and Skipwith 
Pon d [6ij ] (lake or pond); Canton Creek [70| and Stony Stream [7(| (stream); Chesapeake 
Bay [71], St. Mark's Estuary [72], Ythan Estuary [73f . and Ythan Estuary with parasites [73 ], 
(estuaries); Benguela [75|, Caribbean Reef - small [7o|. and North-eastern U.S. continental 
shelf [77j (marine). 

As discussed above, the food webs produced by the evolutionary process in this model 
are relatively small, with average diversity S ~ 15 for the QSS core communities and 
approximately 50 for the full communities. As seen in Fig. [T31 the directed connectance C 
is somewhat smaller (C ~ 0.08 for core and C ~ 0.06 for full communities) than the input 
connectance, c = 0.1. However, this is partly because it is calculated in the conventional 
way as L/S 2 [H3], rather than as L/S{S — 1). The average linkage density is also quite small 
(z = 1.17 for core and z = 2.58 for full communities) compared to most of the real food 
webs that have been documented. For comparison, the seventeen documented food webs 
have S ~ 69 with a range from 25 to 155, C ~ 0.13 with a range from 0.03 to 0.32, and 
z ~ 6.9 with a ran ge f rom 1.6 to 17.8. However, Williams and Martinez' niche model for 
food-web structure [6l( leads to a scaling hypothesis for the degree distributions (exact for 



the niche-model prey distribution and approximate for the predator distribution) [78|, [79|, [80j : 

2zp(k)=p k (k/2z) , (15) 

where k can be either the generality, the vulnerability, or the total degree (with different 
scaling functions p for each). As a consequence, cumulative degree distributions can also 
be rescaled by simply dividing the argument by 2z. Using this niche-model result as a 
general scaling hypothesis [52j, we can compare the degree distributions for the communities 
generated by our model with scaled data for the seventeen real food webs. 
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FIG. 14: (Color online.) (a) Unsealed degree distributions for the QSS core communities, showing 
number of prey (solid), number of predators (dashed), and total degree (dot-dashed). The distri- 
butions were individually normalized for each community, and then averaged over all communities. 
The error bars are standard errors, based on the spread between the twelve simulation runs. Aver- 
age linkage density is z = 1.17. The inset on log-linear scale shows that these distributions decay at 
least exponentially with increasing argument, (b) Same as (a), for the full communities. Average 
linkage density is z = 2.58. (c) Histograms over individual sampled communities of the correlation 
coefficient between a species' numbers of predators and prey. Black filled with light gray (yellow 
online): QSS core communities. Medium gray (red online): full, connected communities. Dark 
gray (blue online) with drop lines: empirical data for seventeen natural food webs. 



The degree distributions for the model, individually normalized for each community, and 
then averaged over all communities, are shown in unsealed form in Fig. fl~47 a) and (b). As 
seen from the insets, all three distributions (prey, predators, and total degree) decay at least 
exponentially for large argument. This property is shared by most real food webs [!59l and 
models [52], [ 
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and indicates that food webs in general are not scale-free networks. 
The behavior for k < 2z is more model and system dependent. 

Histograms for the correlations between a species' generality and vulnerability are shown 
in Fig. fT4Tc). together with data for the seventeen real food webs. The average correlation 
coefficients are —0.53 for core communities and —0.40 for full communities, more negative 
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FIG. 15: (Color online.) Scaled, cumulative degree distributions for the model, compared with 
results for the seventeen empirical food webs. Data for individual, real food webs are shown in 
the background as isolated data points. Average data for the simulated and real webs are shown 
as data points with error bars indicating standard error, connected by heavy lines. Dark gray 
with circles (blue online): simulated core communities. Medium gray with squares (red online): 
simulated full communities. Black with diamonds: empirical communities. Details of the averaging 
are discussed in the text, (a) Number of predators (vulnerability, or outdegree). (b) Number of 
prey (generality, or indegree). (c) Number of predators plus prey (total degree). 

than, but not inconsistent with, the real food-web average of —0.23. 

Due to the relatively small size of the empirical food webs, it is difficult to extract infor- 
mation from scaled probability densities, which will include empty bins unless a very large 
bin size is chosen. This problem is avoided by instead studying the cumulative distribution, 
P(x) = Prob [k/(2z) > x] [H3|. These are shown in Fig. [TS] for the seventeen empirical webs, 
together with averaged results for the empirical webs and the model. The averaged curves 
were generated by first scaling the degrees (predators, prey, and total) for each web by the 
value of 2z for that particular web, binning the results in bins of width 0.2, and normal- 
izing the binned histogram by the species richness S of that same web. The individually 
scaled, binned, and normalized histograms were then averaged over all webs (seventeen for 
the empirical webs and 131 072 x 12 for the model webs) and finally integrated to produce 
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the average scaled cumulative distributions. For the model, results are shown both for the 
full and core communities. 

Despite their smaller average diversity and linkage density, the scaled cumulative degree 
distributions for both types of model communities, which are shown in Fig. [Jj)J f an weu 
within the range of the empirical data. For scaled degrees below about 1.5, the averaged 
distributions for the core communities are in quite reasonable agreement with those for the 
empirical food webs. The largest deviations are for total degree (Fig. fToT c)) at small values of 
the scaled argument. This is possibly due to the significantly stronger negative correlations 
between the numbers of predators and prey (vulnerability and generality) for species in the 
model communities, compared to the empirical ones. (See Fig. IT4Tc).) For scaled degrees 
above approximately 1.5, the averages for the empirical webs are significantly above those 
for the model webs and appear approximately exponential in the tail (see insets in Fig. [T5l) . 
We note, however, that this behavior in the averages is caused by a small number of webs; 
the remaining empirical webs have no species of such highly above-average degrees. Overall, 
the scaled, cumulative distributions for the present model agree with the empirical data to 



about the same degree as the niche model |78j and the speciation model of Rossberg et 
al. [52|, 82|. However, it must be admitted that much potentially valuable detail about the 
food-web structure is lost in calculating the cumulative distributions. 

A measure of the hierarchical structure of a food web is given by the proportions of basal 
species (species with no prey, supported only by the external resource), intermediate species 
(that have both prey and predators) and top species (that have no predators). As seen 
from the histograms in Fig. [THJ there are wide variations from community to community 
for each class of species. The high proportion I of intermediate species seen for the full 
communities in Fig. [TBTb) (/ = 0.76) is consistent with most of the real food webs analyzed, 
including Caribbean Reef (I = 0.94), U.S. Shelf (I = 0.94), Benguela (I = 0.93), Scotch 
Broom (J = 0.92), Skipwith Pond (I = 0.92), Coachella Valley (I = 0.90), Little Rock Lake 
(I = 0.86), El Verde Rainforest (I = 0.69), St. Marks Seagrass (I = 0.69), St. Martin Island 
(/ = 0.69), and Bridge Brook Lake (/ = 0.68), as well as simulations of the Web World 
model [6(j and mean-field analysis of a Lotka-Volterra predator-prey model with evolution 
and competition [83[ ■ The core communities, on the other hand, have their member species 
much more evenly distributed among the three classes, as shown in Fig. [TBTa). Real food 
webs with more even distributions are Canton Creek (B = 0.53, / = 0.22, T = 0.25), 
Stony Stream (B = 0.56, / = 0.27, T = 0.17), and Chesapeake Bay (B = 0.16, / = 0.52, 
T = 0.32). Thus, at least in this model, and somewhat counterintuitively, on average the 
intermediate species are both the least stable and the most numerous. 

We conclude these comparisons of the structures of food webs evolving in the present 
model to data for real food webs with some words of caution. The model was primarily 
developed to study the long-time temporal fluctuations in diversity and population size, and 
it was in no way tuned to produce realistic community structures except insofar as they are 
generally food-web like. This community structure does, however have profound effects on 
the dynamics, as it is responsible for the difference between the power-law exponents for the 
lifetimes of species and communities in this predator-prey model. One should also note that, 
while the lack of correlations between the properties of closely related species was shown 
to have little effect on the long-time dynamics j2lj . we have not considered whether such 
correlations might affect the detailed community structure. Most likely they would, at least 
to some degree. In addition, the growth of real food webs involve a number of mechanisms, 
besides just evolution and local population dynamics. Most notably, immigration plays 
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FIG. 16: (Color online.) Histograms for the proportions of basal, intermediate, and top species. 
The vertical, dashed lines in matching grayscale (color) represent the mean proportion for each 
class, (a) Core communities. Mean proportions: 0.33 for basal, 0.31 for intermediate, and 0.38 for 
top species, (b) Full communities. Mean proportions: 0.14 for basal, 0.76 for intermediate, and 
0.13 for top species. 

an important role in the occurrence of new members of a spatially localized community. 
However, due to the large differences between parent species and mutants in the current 
model, our evolution mechanism may also been seen as covering immigration. 

We also note that the simulated communities have significantly smaller diversities, con- 
nectances, linkage densities, and population sizes than the real communities, so that our 
comparisons had to rely heavily on scaling arguments. It would therefore be desirable in 
future studies to increase the connectance, number of potential species, and resource size to 
bring the results closer to the realistic range. 



V. SUMMARY AND CONCLUSIONS 



In this paper we have studied in detail an individual-based predator-prey model of bio- 



logical coevolution, based on the simplified version of the tangled-nature model [16|, [T7J, [1 
that was introduced in Ref. [HI]. Selection is provided by a population-dynamics model in 
which the reproduction probability of an individual of a particular species depends nonlin- 
early on the amount of external resources and on the population densities of all other species 
resident in the community. New species appear in the community through point mutations 
in a genome consisting of a string of L bits. 

In the mutation- free limit, the mean fixed-point population sizes and stability properties 
of any A/"-species community can be obtained exactly by linear stability analysis. While 
the universal competition effect that enables this analytical treatment is not very realistic, 
the exact solutions make the model ideal as a benchmark for more realistic, but also more 
complicated, models. A preliminary discussion of two more realistic models is found in 
Ref. Q. 

In the simulations presented here, we used L = 20 for a total of 2 20 = 1 048 576 potential 
species. In order to study the statistically stationary properties of the model, we performed 
long kinetic Monte Carlo simulations over 2 25 = 33 554 432 generations. By studying the 
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stationary fluctuations we hope in the future to gain an understanding of the the system's 
sensitivity to external perturbations in a way analogous to the fluctuation-dissipation rela- 



tions of equilibrium statistical mechanics [30 



Qualitatively, many of the statistical properties of this model are similar to those of the 



related, mutualistic model studied in Refs. [19|, |20j, |21|. These include approximate 1/ f 
noise in power spectra (PSDs) of diversity and population sizes (f~ a with a ~ 1.29), and 
power-law distributions for the lifetimes of individual species, as well as of the durations of 
evolutionarily quiet periods, corresponding to QSS communities. However, in contrast to the 
mutualistic model, the power-law exponents for the species lifetimes and QSS durations are 
different: t~ T1 with n ~ 2 for the former (consistent with a stochastic branching process j53j) 
and t~ T with r « —1.07 for the latter. In Ref. [28] it was speculated that the exponent values 
a — 1, T\ — 2, and r = 1 are consistent with predictions for a zero- dimensional extremal- 



dynamics model [85|, |86|. However, this speculation is not consistent with the presumably 
more accurate estimate for a presented here, and it seems advisable to be rather skeptical 
about any mapping of the current model onto a simple statistical-mechanical extremal- 
dynamics model. 

It is probably more fruitful to consider why t% and r are different for the current predator- 
prey model, while they coincide for the corresponding mutualistic model. Here we believe the 
answer lies in the different community structures in the two models. While communities in 
the mutualistic model are tightly knit with all positive interactions (see Fig. 10 of Ref. 29J), 
the communities in the present predator-prey model take the form of simple food webs, 
which are much more resilient toward the loss of a single or a few species. In this sense, this 
predator- pre y model is in much better agreement with real food webs than the mutualistic 
model @, m, H [58 . 



Our comparisons of the structure of the communities generated by our model with real 
food webs show both similarities and differences. A difficulty with such comparisons is 
the large differences in diversity and linkage density between the simulated and real webs, 
which necessitated heavy use of scaling arguments in the quantitative comparisons. For 
comparison of detailed community properties, simulations of models with higher diversity 
and connectance are therefore desirable in the future. Nevertheless, our model presents 
a synthesis of long-term evolutionary dynamics and food-web like communities within an 
individual-based framework of integrated population dynamics and evolution, that should 
provide a sound basis for more refined models in the future. 
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FIG. 17: (Color online.) Normalized correlation functions for the matrix elements Mjj for the 
two schemes described in the text, each based on a single realization of M for L = 6, 8, and 9. (a) 
The scheme described in Ref. jl7T |. (b) The modified scheme introduced here. 

APPENDIX A: MATRIX ELEMENTS FOR LARGER GENOMES 

Here we describe an improved version of the method introduced by Hall et al. 16|, 17] to 
produce pseudorandom matrix elements Mjj for values of L that are too large for the full 
2 L x 2 L matrix M to fit into computer memory. This method permits the matrix elements 
for a given community to be generated and retained only as needed. We first present Hall 
et al.'s method, point out some problematic features, and then present our modifications. 

Let S(7) be the string of binary digits corresponding to the decimal species label /. This 
bit string has length L, so there are 2 L different strings. To generate the matrix element 
Mu, one first generates a new string of the same length, S(J, J) = S(/)XORS(J), where 
XOR is the logical exclusive or operator. From this bit string is generated the corresponding 
new decimal index K (S(J, J)). Next one creates two one- dimensional arrays, X and Y, each 
of 2 L random numbers between —1 and +1. (For simplicity let the starting index for the 
arrays be zero.) Since S(J, J) is symmetric in I and J, asymmetric pseudorandom matrix 
elements are generated as 

M u = [X(AT(S(/,J)))+Y(J)]/2. (Al) 

(Hall et al. instead use the product of the two random numbers, which gives a pseudorandom 
number with a slightly different distribution.) 

The problem with this method is that it produces strong correlations along the columns 
of M since the second of the two random numbers that produce Mu is the same for all 
I at the same J. In Fig. [TTT a) we show the resulting correlation function for this scheme 
as a function of the Hamming distance between the pairs of bit strings involved in two 
different matrix elements. Regardless of L, significant correlations are seen for Hamming 
distances less than five. A discussion of how to calculate such correlation functions is found 



in Ref. 21 



To reduce these correlations between matrix elements involving closely related genotypes, 
we here modify the scheme as follows. We extend array Y to a length of 3 x 2 L and define 
M u as 

Mu = [X (K (S(J, J))) + Y(K (S(J, J)) + 2( J + 1))] /2 . (A2) 
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FIG. 18: (Color online.) Histogram for the pseudorandom Mjj produced by the scheme introduced 
here for L = 9, corresponding to 262 144 individual matrix elements. The straight, black lines 
represent the triangular shape of the theoretical distribution. 

The addition of K (S(I, J)) in the index of Y ensures that the matrix element depends in an 
erratic fashion on both / and J, while the term linear in J ensures that M is not symmetric. 
The resulting correlation function is shown in Fig. fTTT b). The correlations for elements 
involving closely related genotypes are strongly suppressed. The correlations for a Hamming 
distance of 2L, which are present in both schemes, are of little practical significance. They 
are caused by the fact that the XOR operation is invariant under simultaneous bit reversal in 
both its arguments and could be removed by adding a linear function of I in the argument 
of X. The probability density for M/j is triangular as expected from simple analytical 
arguments, and the whole interval from —1 to +1 is well sampled, even for relatively small 
L. This is shown in Fig. [T8l 
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